#! usr/bin/perl -w
use strict;

# print results on two output files
sub output
{
	my($ref_list_in,$ref_list_out)=@_;
	my @list_in=@$ref_list_in;
	my @list_out=@$ref_list_out;

	my @uniq_list_in = keys %{{ map { $_ => 1 } @list_in}};	
	my @uniq_list_out = keys %{{ map { $_ => 1 } @list_out}};

	open(OUTPUT1,">SNP_in.txt") or die "cannot create outputfile";
	open(OUTPUT2,">SNP_out.txt") or die "cannot create outputfile";

	foreach my $item1 (@uniq_list_in)
	{
		print OUTPUT1 "$item1\n";
	}		
	
	foreach my $item2 (@uniq_list_out)
        {
                print OUTPUT2 "$item2\n";
        }

	close OUTPUT1;
	close OUTPUT2;

	print "\nsee independent SNPs in the output file \"SNP_in.txt\" and others in \"SNP_out.txt\" \n\n";

}

# return a boolean (if an item is in the table or not)
sub inTab
{
	my($item,$ref_tab)=@_;
	$item=shift;
	my @tab=@$ref_tab;

	my @in=grep(/^$item$/ig,@tab);

	if(@in eq 0)
	{
		return 0;
	}
	else
	{ 
		return 1;
	}
}

# return position of an item
sub position
{
	my($item,$ref_tab)=@_;
        $item=shift;
        my @tab=@$ref_tab;

        my $pos="";

        for (my $i=0; $i lt (scalar@tab); $i++)
        {
                if($tab[$i] eq $item)
                {
                        $pos = $i;
                }
        }

        return $pos;	
}

# delete an item from the table
sub delTab
{
	my($item,$ref_OldTab)=@_;
	my @OldTab=@$ref_OldTab;
	
	my $pos=&position($item,\@OldTab);
	my @NewTab=();

	for (my $i=0; $i lt (scalar@OldTab); $i++)
	{
		push(@NewTab, $OldTab[$i]) if ($i ne $pos);
	}	

	return(@NewTab);
}

# return two list : list_in and list_out(SNPs'list to delete)
sub selection
{
	my($ref_list_in,$ref_data,$ref_ToListOut)=@_;
	my @list_in = @$ref_list_in;
	my %data = %$ref_data;
	my @ToListOut=@$ref_ToListOut;

	my @list_out=();

	while (my($key,$val) = each(%data))
	{	
		#@tab >=1 if exists in list_in
                my @key_in=grep(/^$key$/ig,@list_in);#search key in list_in (list of independent SNP)
                my @val_in=grep(/^$val$/ig,@list_in);#search val in list_in

                #@tab >=1 if exists in list_out
                my @key_out=grep(/^$key$/ig,@list_out);#search key in list_out
                my @val_out=grep(/^$val$/ig,@list_out);#search val in list_out

		if((@key_in eq 0) and (@key_out eq 0) and (@val_in eq 0) and (@val_out eq 0))#key and val aren't on a list
		{
			push(@list_in,$key);
			push(@list_out,$val);
		}
		elsif(@key_in eq 1)#key already in list_in  
		{
			push(@list_out,$val) if((&inTab(&val,\@list_out)) eq 0);
		}
            	elsif(@key_out eq 1)#key already in list_out
		{
                        push(@list_out,$val) if((&inTab($val,\@list_out)) eq 0);
			@list_in=&delTab($val,\@list_in);
		}
		elsif(@val_in eq 1)#val already in list_in
		{
			push(@list_out,$key) if((&inTab($key,\@list_out)) eq 0);	
		}
		elsif(@val_out eq 1)#val already in list_out
		{
			push(@list_out,$key) if((&inTab($key,\@list_out)) eq 0);
			@list_in=&delTab($key,\@list_in);
		}
		#print "Key $key has value $val\n";
	}
	
	# add SNPs to list_out
	foreach my $item (@ToListOut)
	{
		if((grep(/^$item$/ig,@list_out)) eq 0)
		{
			push(@list_out,$item);
		}
	}
 
	return(\@list_in,\@list_out);
}

sub load_data
{
	# parameters
	my ($filename,$t)=@_;

	my $firstLine=1;
	my @list_in=();
	my @ToListOut=();
	my %data=();

	open(LDFILE, $filename) or die "cannot open $filename";
	while(my $line=<LDFILE>)
	{
		chomp($line);
		$line=~ s/ +/ /gi;#separate rows with a whitespace only
		my @fields=split(/ /,$line);

		if($firstLine!=1)#not header
		{		
			# collect fields of interest
	                my $R2=$fields[(scalar@fields)-1];# field R2
        	        my $SNPid1=$fields[3];# SNP1
                	my $SNPid2=$fields[6];# SNP2
			
			# r2 < threshold : SNP independent so put in list_in
			if ($R2<$t)
			{
				my @match1=grep(/^$SNPid1$/ig,@list_in);#search SNPid1 in list_in
                                my @match2=grep(/^$SNPid2$/ig,@list_in);#search SNPid2 in list_in
	
				push(@list_in,$SNPid1) if(@match1 eq 0);
				push(@list_in,$SNPid2) if(@match2 eq 0);			
			}
			# else r2 > threshold : put in table to analyse
			else
			{
				if(!exists($data{$SNPid1})) 
				{$data{$SNPid1}=$SNPid2;}
				else
				{push(@ToListOut,$SNPid2);}
			}
		}

		$firstLine=0;

	}
	close LDFILE;

	# SNP selection
        my($ref_new_list_in,$ref_list_out)=&selection(\@list_in,\%data,\@ToListOut);	
	my @new_list_in=@$ref_new_list_in;
	my @new_list_out=@$ref_list_out;
	
	return(\@new_list_in,\@new_list_out)
}


sub main{

if (scalar@ARGV != 2)
{
	print "\nExpected 2 parameters. \$ perl snpIndep.pl <LDfilename> <threshold>\n";
}
else
{	
	# parameters
	my $LDfilename=$ARGV[0];
	my $threshold=$ARGV[1];
	
	# load data and SNP selection
	my($ref_list_in,$ref_list_out)=&load_data($LDfilename,$threshold);
	my @list_in=@$ref_list_in;
	my @list_out=@$ref_list_out;

	# output
	&output(\@list_in,\@list_out);
}


exit;
}

main();
